Jamming and pattern formation in models of segregation 
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We investigate the Schelling model of social segregation, formulated as an intrinsically non- 
equilibrium system, in which the agents occupy districts (or patches) rather than sites on a grid. We 
show that this allows the equations governing the dynamical behaviour of the model to be derived. 
Analysis of these equations reveals a jamming transition in the regime of low- vacancy density, and 
inclusion of a spatial dimension in the model leads to a pattern forming instability. Both of these 
phenomena exhibit unusual characteristics which may be studied through our approach. 
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I. INTRODUCTION 
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Forty years ago Thomas Schelling proposed a simple 
model of social segregation in which agents of two differ- 
ent types are placed at random on a grid, before being 
allowed to move according to a desire to be close to other 
agents of the same type [l( . This seminal work played an 
important part in the development of the new scientific 
field of social simulation, in which sociological problems 
are studied through computational models. 

To the physics community, there is an immediate sim- 
ilarity between the rules of the Schelling model and some 
of the simple dynamical models of statistical physics and 
the tools of this field have been used to provide quantita- 
tive insight into the behavior of the model [2| . There have 
also been several attempts to make direct links to equi- 
librium spin models though care must be taken as 
agents in the Schelling model are subject to kinetic con- 
straints which mean that the dynamics cannot be viewed 
as a simple energy minimisation process [2|. In what 
follows we will show that certain analytically tractable 
Schelling-like models of segregation exhibit a range of in- 
teresting physical phenomena. 

In Schelling's original work and most subsequent 
studies , the city in which the agents reside is modeled 
as a 2-dimensional grid. We propose instead to study a 
model in which the basic object of interest is a district 
(patch) containing multiple residences. This modifica- 
tion improves the social realism of the model - a city 
is much better described as a collection of districts and 
suburbs with their own ethnic character than as a sim- 
ple grid. Moreover, patch models are naturally amenable 
to analytical approaches, for example in ecology Q. In- 
deed, a patch variation of a Schelling-like model has been 
considered before [4], though that work has a very dif- 
ferent flavor to our own, being primarily concerned with 
applying the techniques equilibrium statistical mechan- 
ics. We choose the dynamical rules of our model to be 
close to those of the lattice-based model of Gauvin et 
al [5|, simulations of which display interesting physical 
behavior that we intend to study theoretically. 

Our analysis is divided into two parts, investigating 
different implementations of the patch-based Schelling 



model. In both cases we consider a large city, divided 
into N patches each containing K residences. For the 
first model we investigate, model A, the patches are rela- 
tively small, but there are very many of them. Through 
enumerating the possible interactions between patches, 
we obtain a description of the model in terms of a de- 
terministic dynamical system. This framework is used to 
investigate a jamming transition present in the model, in 
which large numbers of agents remain stuck in unfavor- 
able states. We then go on to analyse a model with a 
spatial structure, model B, and consider the alternative 
limit of very large patches. The behavior in this case is 
rather different, with the model exhibiting pattern for- 
mation driven by antidiffusion. 



II. MODEL A 

We begin by considering the situation where K is rela- 
tively small; each patch represents a local neighborhood 
containing only a few residences. Initially, the city is 
randomly populated with equal numbers of agents of two 
different types, which we call A and B, with a fraction 
p of the residences left vacant. At each time step two 
residences are chosen at random from the whole city. If 
the first contains an agent and the second is vacant, then 
that agent is given the opportunity to move to the vacant 
residence. They take up this offer only if the number of 
agents of the opposite type in the destination patch is 
less than a threshold T. 

The contents of patch i at time t is encoded in the num- 
bers dj(i), bi(t) and Vi(t) of A agents, B agents and va- 
cancies it contains. The state of the system as a whole is 
then specified (up to trivial reordering of the residences) 
by the quantities 
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Fa,b,v(t) — — / ] ^g, ai (t)^b,bi{t)^v,Vi(t) > 
i=l 



(1) 



giving the fraction of patches in state (a, b, v) at time t. 
Our theoretical work is based on an analysis of the time 
evolution of these quantities when the number of patches 
is very large. 
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FIG. 1: Equilibrium interface density Xoo as a function of 
vacancy density p, for patches of size K = 9 and tolerance 
parameter T = 3. Theoretical prediction (solid line) from 
equation ([2]) is compared to results from stochastic simula- 
tions on N = 10 3 patches (black circles) and on the lattice 
(grey squares), averaged over 100 samples. Error bars are one 
standard deviation. 



The first step is to consider the possible changes to a 
given F a ^, v which can occur in one time step. Suppose, 
for example, that at time t an A agent in a patch with 
state (a, b, v) is selected to move to a vacancy in a patch 
with state (a',b',v'). This event occurs with probability 

F a ,bAt)^F a ,, b , iV ,{t)^e(T - b') . 

The factors in this expression are explained as (i) the 
probability of choosing a patch in state (a, b, v), (ii) the 
probability of selecting an A agent from that patch, (iii) 
the probability of choosing a patch in state (a 1 , b 1 , v') for 
the destination, (iv) the probability of selecting a vacant 
residence in the destination patch, and lastly (v) a Heav- 
iside 9 function imposing the requirement that the des- 
tination patch contains fewer than T agents of type B. 
As a result of this interaction, the values of F a ^ v (t) and 
Fa',b',v'(t) would decrease by 1/N, whilst F a -i,b,v+i(t) 
and F a / + x,b,v'-i(t) would increase by 1/N. The effects of 
moving a B agent are computed in the same way. 

In the present model, the patches are well-mixed in 
the sense that each patch is equally likely to interact 
with every other; it follows that when the number of 
patches is very large, it is sufficient to consider only the 
average over all possible interactions. Rescaling time by 
a factor of 1/N, we take the thermodynamic limit N — > 
oo, in which the random quantities F a ^_ v (t) may be well 
approximated by continuous deterministic functions of 
rescaled time. In this formalism the possible changes that 
can occur in one time step discussed above are translated 
into a deterministic differential equation which exactly 
describes the behavior of the system in the limit N — > oo. 
Summing over the possible interactions gives 
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where the contributions from each type of interaction are 



Ra = 



Ra = 



K 

y_ 

K 
K 
K 



w-1 



K 


b) 


K 


a) 


F a ,b,v(t) 


a " 

K- 


F a ,b,v(t) 


b " 
K- 



The £ ± used here are step operators which alter the func- 
tions they act on through the addition or subtraction 
of 1 to their argument; for example, £~£^[F a ^ 7V (t)] = 

To monitor the emergence of segregation in the model, 
we measure the fraction of pairs of neighbouring agents 
of different types, a statistic commonly referred to as the 
interface density 0, [H, H| . A patch in state (a, b, v) has 
(a + b)(a+b— 1)/2 distinct pairs of agents, of which ab are 
(a, b) pairs. The interface density x is found by summing 
over all patches, or alternatively by the formula 



5(t) - J2 *"«.*.«(*) 



2ab 



a,b,v 



(a + b)(a + b-l) 



In simulations of lattice-based versions of the model 
(without patches), behavior suggesting a phase transi- 
tion in interface density has been observed as the vacancy 
density p is lowered; below a critical value, there are not 
enough vacancies to facilitate the movement of agents 
to a segregated state and the system appears 'jammed' 

Q . We can confirm the existence of this transition in 
the patch model under investigation here by numerically 
integrating @ to find the final state of the system. 

Starting from a well-mixed initial condition for the 
F a ,b,v (chosen to be equivalent to the large N limit of a 
random initial configuration of agents in the microscopic 
model) we numerically integrate ^ using the forward- 
Euler scheme. Figure 1 shows the equilibrium value Xoo 
of interface density as a function of p for patches of size 
K = 9 and a tolerance parameter T = 3. 

For low values of vacancy density the segregated steady 
state becomes inaccessible to the dynamics started from 
a well-mixed initial condition, and the system finds a dif- 
ferent (well-mixed) equilibrium; we estimate the critical 
point for the transition to be p c « 0.079615. Further 
information about the system can be gained through a 
linear stability analysis of the steady states reached from 
both well-mixed and segregated initial conditions. In the 
jammed regime we find two stable equilibria, one well- 
mixed and one segregated, moreover, the stability of both 
states increases with p. Past the transition point only a 
single, segregated, steady state can be found. 

In figure 1, the analytical result is compared with re- 
sults from simulations of the patch model with N = 10 3 
patches, averaged over 100 samples. Also shown are the 
results from simulations of a lattice based version of the 
model [5j, where the jamming transition occurs at a dif- 
ferent point. We should point out that the transition is 



3 




600 



t 



FIG. 2: Comparison between the deterministic theory (black 
line) and the results of ten simulation runs (gray lines) for va- 
cancy density p = 0.08, just above the point of the jamming 
transition. Patches of size K = 9 and a tolerance parame- 
ter of T — 3 were used for both data sets. The stochastic 
simulations were performed on a system of N = 10 3 patches. 



not unique to the values of K and T we have chosen; fur- 
ther numerical results suggest that the conditions K > 3 
and T < K/2 are sufficient. 

Beyond confirming the existence of the transition, we 
are also able to probe the behavior of the system near 
the critical point. Figure [5] shows the deterministic dy- 
namics for a value of vacancy density just above critical, 
p = 0.08. Also shown in that figure are several sample 
results from stochastic simulations of a system of size 
N = 10 3 . The data from simulations show very little 
noise, although both the final outcome (either jammed 
or unjammed), and the moment in time that unjamming 
takes place, appears random. Ordinarily, one might ex- 
pect that the persistence of a metastable state followed 
by sudden relaxation to equilibrium is a stochastic effect 
which would not be captured by a naive deterministic 
theory. In the present case, however, we see that the de- 
terministic theory does indeed display the same behavior, 
with the effect of stochasticity and finite size mainly lim- 
ited to fluctuations in timing. These effects are reduced 
in larger system sizes. 

There is an analogy [!, Q between the lattice based 
version of the segregation model and kinetically con- 
strained models which are used as simplified proxies in 
the study of glasses and granular media [Tof. Kineti- 
cally constrained models have been intensively studied 
over the last decade and they exhibit a rich phenomenol- 
ogy, including a jamming transition which is both 
discontinuous in its order parameter and features expo- 
nentially diverging relaxation time. 

In the present model it is also the case that, as vacancy 
density is further lowered towards the critical value, the 
waiting time until the system is freed from the metastable 
jammed state increases. In analogy with critical slowing 
of magnetization in the Ising model [Til ], we introduce 
the following relaxation time for interface density: 
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dt. 



(3) 




viewed as a function of p — p c . As shown in Figure [3J this 
quantity does not grow exponentially as in some kineti- 



FIG. 3: Log-log plot of the duration of the metastable state, 
measured by R given in equation ([3]), against p — p c . Circles 
show the results of numerical integration of p]l. the dashed 
line is a power-law with exponent —1/2. 



cally constrained models, but rather exhibits power-law 
behavior with exponent —1/2. 



III. MODEL B 

We now introduce a simple modification to the above 
model which will allow us to develop a theoretical de- 
scription of the emergence of spatial patterns of agents. 
We take the N patches considered earlier and arrange 
them in a 1-dimensional lattice; from now on, agents are 
restricted to move only to patches which neighbor them 
in the lattice. The other rules of the model remain the 
same, though to aid the analysis we smooth the step func- 
tion, introducing Q K (x) = (l+tanh(2x/«;))/2, for a small 
parameter k. 

For this spatial version of the model we employ a dif- 
ferent analysis in which we keep the number of patches 
finite and take the limit of large patch size K — > oo. 
In a sociological context this model represents a city di- 
vided into several districts, each still containing a large 
number of residences. A similar model featuring agents 
with no preferences was considered in [15| and the same 
techniques apply here with the simple addition of the K 
threshold term. 

The state of the system at a given time is specified by 
the vectors a = (oi, . . . a/v) and b = [b±, . . . 6jy) giving 
the numbers of A and B agents in each patch. Note that 
there is no need to keep track of vacancies as a% + bi + Vi = 
K for all i. The dynamics of the model are determined 
by the transition probabilities P{a, b\a', b'), giving the 
likelihood of moving from state (a',b f ) to state (a, 6) in 
one time step. Changes to the system result from the 
movement of agents between neighboring patches; the 



4 




transition probabilities for A and B agents have the forms 

a,; Vi 



t = 5 



P(ai - l,aj + l\ai,a,j) = 



e K (6, - t) : 



P(6 i -l ! 6 J + l|6 i ,6 i ) =*|i_ i | 1 



NK2K V 7 ' 



(4) 



where i and j are neighbors in the lattice, and we have 
listed as arguments only those entries of a and b which 
change. Writing 7r t (a,b) for the probability of finding 
the system in state (a, 6) at time t, we have the equation 



7r t+ i(o,6) -7r f (o,6) 

= (P(a,b\a' ,b')ir t (a' ,b') - P(a' ,b'\a,b)n t (a,b) 



Multiplying by a, and summing over all states, we obtain 
a difference equation for the average number of A agents 
in patch i: 

{a l ) t +i - (ai)t = 

^P(a^ + 1, aj — l\cti, Oj) — P(ai — 1, Oj + l|aj, a 

where the notation j £ i indicates that the sum is over 
all patches j which are nearest neighbors of patch i. A 
similar expression exists for B agents. We now rescale 
time by a factor of 1 / K and introduce 



K 



and t = — . 

K 



Taking the limit K — > oo transforms the difference equa- 
tions for (aj)t and (bi)t into a pair of coupled reaction- 
diffusion equations: 



oti = 1&k(t - ft)Aai - ajA[7,-6 K (r - ft)] 
ft = TiQ^-aOAft-ftA^e^r-a;)] 



(5) 



where A is the discrete Laplacian. 

By taking the limit of a large number of agents per 
patch, the emergence of segregation in this model is re- 
duced to the question of stability of the homogeneous 
state oti = Pi = (1 — p)/2 for all i. Diagonalizing the 
Jacobian of ([5]) at this point, we find that there is once 
again there is a jamming transition: for fixed t < 1/2 
and the homogeneous state is stable for small p, 

becoming unstable as p approaches 1 — 2r. 

In numerical simulations starting with a uniform dis- 
tribution of agents in the un-jammed phase, the model 
rapidly forms a distinct pattern of alternate patches filled 
by agents of different types. The emergence of these pat- 
terns can be studied theoretically by setting p = 1 — 2r 
and analysing the stability of a homogeneous distribution 
of agents (that is, a« = ft = (1 — p)/2 for all i). 

In the analysis of typical pattern forming systems, it is 
normal to study the continuum version of the reaction- 
diffusion equations [l6l |. The homogeneous state may be 




FIG. 4: (Color online) Upper - two snapshots taken from 
a simulation of the spatial patch model on a 1-dimensional 
toroidal lattice of N — 10 patches, with patch size K — 10 6 , 
vacancy density p — 0.1 and step parameter k = 0.1. The 
fraction of agents of type A and B in each patch are shown 
by darker (red) and lighter (blue) areas, respectively, with va- 
cancy density in white. Lower - time evolution of the average 
absolute difference in agent densities (10!) = jf ^ \cn — Pi\ 
(solid line), alongside the linearised result (dashed line) from 
the theory. 



found to be unstable with respect to perturbations which 
are periodic in space; the eventual shape of the pattern 
which forms is then given by the most unstable wave- 
length. In the present model we find a different behavior, 
with instability on the scale of lattice spacing making it 
inappropriate to take the continuum limit. 

The linearised form of ([5]) can be decoupled by con- 
sidering vacancy densities 7, and a conjugate variable 
Q = OLi — ft. Near the homogeneous fixed point we find 



pO--pY 



c 



(6) 



AO 



From the first line above we infer that, near the ho- 
mogeneous state, the vacancies will exhibit diffusion, in- 
dependent of the behavior of the agents. For the con- 
jugate variables Ci the behavior is opposite: note that 
the coefficient in the second line is typically large and 
negative, indicating rapid anti-diffusion. By taking an 
initial condition in the form of an alternating perturba- 
tion on = (l-/9)/2 + e(-l) l andft = (l- ( o)/2+e(-l) l+1 
for a small positive e, we find that the lattice dependence 
drops out, giving 

li {t) = p, and Ci(t)=2e(-l) i e-(P- ei ^ £l ) t . 

Together these facts give a description of the emer- 
gence of an alternating pattern of patches dominated by 
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agents of different types, which develop without altering 
the distribution of vacancies. Initially the pattern will 
grow exponentially quickly before being limited by the 
effect of patches becoming saturated with one type of 
agent, at which time the stochastic model will deviate 
from the results of the anti-diffusive linearised system. 
This picture is confirmed by comparison to numerical 
simulations, as shown in Figure UJ 

IV. DISCUSSION 

The models introduced and analysed above are sim- 
ple coarse-grained modifications of the well-known lat- 
tice based segregation models introduced by Schelling. 
We have demonstrated that the well-mixed patch model 
reproduces the jamming transition observed in simula- 
tions of lattice models, and moreover that it is described 
by deterministic equations which are exact in the ther- 
modynamic limit. One might expect that the unusual 
character of the jamming transition is intrinsically linked 
to the effects of stochasticity or spatial dimension. Our 
analysis shows that this is not the case; the behavior is 
captured by the deterministic dynamical system @. 

It should be noted that, despite the similarities in 
the nature of the models and of the characteristics of 
the transition, the jamming found here is possibly a dif- 
ferent phenomenon to that occurring in kinetically con- 
strained models. Specifically, in those models jamming 
often refers to the (probabilistic) existence of an exten- 



sive fraction of frozen spins |llTll3| . The analysis we have 
presented links the occurrence of jamming in a stochastic 
model to a dynamical transition in a particular determin- 
istic system. As such, we cannot draw conclusions about 
the behavior of individual agents in the thermodynamic 
limit. Another possible analogy is with quasi-statio nary 
states in models with long range interactions [13, fl8| . 
Whatever its precise nature though, the jamming tran- 
sitions in models of segregation are certainly interesting 
phenomena which are worthy of further investigation. 

Introducing a spatial dimension to the model and 
studying the limit of large patch size, we observe pat- 
tern formation driven by anti-diffusion. Once again, the 
model is simple enough that the theoretical analysis pro- 
vides a very complete picture of the mechanisms behind 
the emergence of segregation. 

These results exemplify the possibilities for theoretical 
physics analyses of social models. Whilst our emphasis 
has mainly been on the novel physics encountered in such 
systems, the predictive power of the theoretical approach 
we use should also be of interest to simulators. 
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